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Abstract. 

Depth distributions of positron-emitting nuclei in PMMA phantoms are 
calculated within a Monte Carlo model for Heavy-Ion Therapy (MCHIT) based 
on the GEANT4 toolkit (version 8.0). The calculated total production rates 
of 1X C, 10 C and ls O nuclei are compared with experimental data and with 
corresponding results of the FLUKA and POSGEN codes. The distributions of e+ 
annihilation points are obtained by simulating radioactive decay of unstable nuclei 
and transporting positrons in surrounding medium. A finite spatial resolution of 
the Positron Emission Tomography (PET) is taken into account in a simplified 
way. Depth distributions of /3 + -activity as seen by a PET scanner are calculated 
and compared to available data for PMMA phantoms. The calculated /3+-activity 
profiles are in good agreement with PET data for proton and 12 C beams at 
energies suitable for particle therapy. The MCHIT capability to predict the 
Z?"*" -activity and dose distributions in tissue-like materials of different chemical 
composition is demonstrated. 
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1. Introduction 

Beams of charged particles, protons and light nuclei, are proved to be very efficient for 
the radiation therapy of deep-seated solid tumours (Castro et al 2004) . Their selective 
and strong impact on a tumour is possible due to a maximum energy deposition at 
the end of their range in matter (the Bragg peak) . According to the statistical data of 
Particle Therapy Cooperative Group (PTCOG 2006) over 48000 patients were treated 
worldwide with proton or ion beams by July 2005. 

Light nuclei, e.g. carbon ions, have an additional potential advantage compared to 
protons associated with their enhanced relative biological effectiveness (RBE) (Kraft 
2000). Treatment schemes using the carbon- ion irradiation demonstrated good clinical 
results (Schulz-Ertner et al 2004, Tsujii et al 2004, Kanai et al 2006). 

In most cases so far (see PTCOG 2006) the dedicated therapy facilities use either 
proton or heavy- ion beams. Since 2001, the Hyogo Ion Beam Medical Center (HIBMC) 



Positron emitting nuclei in proton and heavy-ion therapy 



2 



is the first facility in the world which operates with two types of beams, protons and 
carbon- ions, in a single location (Hishikawa et al 2004, Mayahara et al 2005). Several 
new facilities are planned or under construction in Europe: HIT in Germany (Haberer 
et al 2004, Heeg et al 2004), CNAO in Italy (Amaldi 2004), ETOILE in France 
(Bajard et al 2004) and MedAustron in Austria (Griesmaycr and Auberger 2004). 
These particle therapy centres will also have proton and carbon-ion beams on their 
sites. Each of the two types of beams or their combination will be available for patient 
treatment. 

Successful particle therapy requires thorough treatment planning, which must 
include optimisation between positive and negative effects of the beam. A required 
dose of radiation has to be delivered to a tumour while sparing surrounding healthy 
tissues and organs at risk. The optimisation of doses requires a reliable method for 
calculating the beam energy deposition in highly heterogeneous human tissues. Up 
to now, mostly deterministic methods have been used by different groups for dose 
calculations in hadron therapy, see e.g. Hong et al (1996), Kramer et al (2000), Jakel 
et al (2001). However, nowadays a Monte Carlo simulation of particle transport in 
human body can be used as an alternative method for dose calculations. Increasing 
power of modern computers makes this method less restrictive in terms of the CPU 
time (Rogers 2006). 

In view of the prospective wide use of both proton and carbon-ion therapy 
in several medical centres worldwide, one should think of employing common 
computational tools for proton and ion therapy based on the Monte Carlo approach. 
We believe that the GEANT4 toolkit (Agostinelli et al 2003) is well suited for this 
purpose. This toolkit was created by an international collaboration of physicists and 
programmers for basic research in nuclear and particle physics. This is an open-source 
project based on modern techniques of programming and visualisation. It is capable 
of simulating a wide range of physical processes in extended media, which makes it 
useful for applications in medical physics. 

The GEANT4 object-oriented toolkit is written in C++ and contains classes 
and methods which provide the basic functions for simulations (GEANT4-Documcnts 
2006), e.g. handling setup geometry, tracking primary and secondary particles, 
calculation of their energy loss and energy deposition in matter, run management, 
visualisation and user interface. On this basis, a specific application can be developed 
for each particular task. The results on validation of the GEANT4 toolkit for proton 
therapy (Jiang and Paganetti 2004) and for therapy with electrons and photons 
(Carrier et al 2004, Larsson et al 2005) have been presented recently. A GEANT4- 
based application for emission tomography (GATE) was developed by the Open-GATE 
collaboration (Jan et al 2004) and successfully used in Monte Carlo simulations 
of commercial PET scanners (Schmidtlein et al 2006, Lamare et al 2006). The 
GEANT4 toolkit is also used for simulations in brachytherapy (Enger et al 2006, 
Perez-Calatayud et al 2006). 

In our recent paper (Pshenichnov et al 2005) we have made a first attempt to 
validate the GEANT4 toolkit of version 7.0 (GEANT4-Webpage 2006) for heavy-ion 
therapy simulations. We have developed a GEANT4-based application for Monte 
Carlo simulations of Heavy-Ion Therapy (MCHIT). The depth-dose distributions 
in tissue-like media, calculated within this model, are in a good agreement with 
experimental data. The relative contribution to the total dose due to secondary 
neutrons was quantitatively evaluated for proton and ion beams. 

In the present work we extend the applicability of the MCHIT model to 
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calculations of secondary /3 + -radioactivity induced by proton and heavy-ion beams 
in tissue-like media. Currently the MCHIT application is based on the version 8.0 
(with patch 01) of GEANT4 toolkit. 

2. Dose and range monitoring in particle therapy by positron emission 
tomography 

The main advantage of proton and heavy-ion therapy consists in the possibility to 
deliver a prescribed dose to a tumour volume while sparing surrounding healthy tissues 
and organs at risk. This requires thorough verification of the proton and heavy- 
ion ranges in patient's body. Tissue inhomogencitics and local anatomical changes 
between fractions of particle therapy may lead to some differences between the actual 
dose delivered to a patient and the prescribed dose. 

Two methods for dose verification in proton and heavy-ion therapy were proposed. 
In the first method, the beams of /3 + -radioactive nuclei 19 Ne, 7\/ 2 = 17.22 s, (Tobias 
et al 1977), n C, T 1/2 = 20.39 min, (Urakabe et al 2001) or 10 C, T 1/2 = 19.255 s, 
(Iseki et al 2004) are injected into the patient, fixed in the proper position, just before 
the therapeutic treatment with the beams of stable nuclei, 20 Ne or 12 C, respectively. 
However, the practical application of this method is limited by high costs of the 
radioactive beams. 

More frequently another method is used which utilises the fact that [3 + - 
radioactive nuclei are produced in fragmentation reactions taking place in tissues 
during irradiation with beams of protons or stable nuclei, e.g. 12 C. In proton therapy 
target fragmentation reactions produce 15 0, T X j 2 = 122.24 s, among other fragments. 
As was early realized (Bennett et al 1975, 1978), this can be used for controlling 
the proton beam localisation in tissues. Later, a similar auto-activation technique, 
but involving 10 C and n C nuclei produced in fragmentation of 12 C projectiles, was 
proposed for dose monitoring in carbon-ion therapy (Pawelke et al 1996, 1997). 

Unstable nuclei 10 C, n C and 15 undergo /3+-decay: A(Z,N) -> A(Z -l,N + 
1) + e + + u c . Due to the three-body kinematics, the energy released in the transition 
of a bound proton into a bound neutron is partly carried away also by a neutrino, 
so that emitted positrons have a continuous energy spectrum. The maximum energy 
of positrons emitted e.g. by 11 C nucleus is 960 keV with the average energy of 386 
keV. Such positrons can travel up to 4 mm in the human body before they stop 
and annihilate on surrounding electrons: e + c~ — > 77. Due to multiple scattering 
of positrons, their path in a tissue is far from the straight line, and their average 
penetration depth is shorter than their actual path. This difference is described by 
a detour factor (Fernandez- Varea et al 1996), which is estimated, for example, for 
1 MeV electrons or positrons in water as ~ 0.5. As a result, most of the positrons 
emitted by 11 C annihilate within ~ 2 mm from their emission point (see also Levin 
and Hoffman 1999). 

In the course of multiple scattering the positrons slow down, so that e + e~ 
annihilation take place practically at rest. Therefore, the angle between the momenta 
of the emitted photons is close to 180°. They are registered by detectors outside the 
patient body. Corresponding reconstruction algorithms make possible to obtain the 
spatial distribution of annihilation points, which should be close to the distribution 
of positron-emitting nuclei. It is further expected that the spatial distribution of 
/3 + -activity induced by proton and heavy-ion beams is strongly correlated with the 
corresponding dose distribution. However, there exists no simple way to express the 
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dose distribution in terms of the /3 + -activity distribution. 

In practice, the problem is solved in a few steps. First, the /3 + -activity 
distributions are measured in advance in experiments with phantoms for different 
beam energies and doses. In particular, such measurements have been performed 
for proton beams in phantoms made of polymethylmethacrylate (PMMA), C5H8O2 
or lucite (p = 1.18 g/cm 3 ) by Oelfke et al (1996) and Parodi et al (2002). The 
/3+-activity measurements for 20 Ne beams in PMMA were made by Enghardt et al 
(1992), and later for 12 C beams by Pawelke et al (1997), Ponisch et al (2004) and 
Parodi (2004). 

Second, the models capable of calculating both the dose and /3 + -activity 
distributions are validated with these data. For example, a phenomcnological model 
for the proton transport incorporating data on isotope production was used by Oelfke 
et al (1996) to calculate the dose and /3 + -activity distributions in homogeneous media. 
Parodi and Enghardt (2000), Parodi et al (2002), Ponisch et al (2004) simulated the 
transport of protons and carbon ions in PMMA phantoms with FLUKA and POSGEN 
codes, respectively. In the present paper we propose to use the GEANT4 toolkit for 
calculating the dose and /3 + -activity distributions. 

Third, the dose applied to a patient in a therapeutic treatment can be verified via 
comparison of the /3 + -activity distribution in the patient body with the distribution 
predicted by the model for the same dose. 

3. GEANT4 physics models used in MCHIT 

We use the version 8.0 (with patch 01) of the GEANT4 toolkit (GEANT4-Webpage 
2006) to build a Monte Carlo model for Heavy-Ion Therapy (MCHIT). Currently this 
model is capable of calculating the three-dimensional distributions of dose and /in- 
activity in tissue-like media. In the present study we use homogeneous phantoms 
represented by a water cube or by a cube made of PMMA. Calculations can also 
be made for phantoms representing bone, liver or muscle tissues. In simulations the 
phantom is irradiated by beams of protons or ions with given beam size, emittance, 
angular convergence/divergence and energy distribution. The phantom material and 
size, as well as all beam parameters can be set via user interface commands individually 
for each run. 

According to the GEANT4 concept, the set of physical models which are relevant 
for a particular problem should be activated by the user via a set of commands. 
A detailed description of physical models included in GEANT4 is given in the 
Physics Reference Manual (GEANT4-Documents 2006). In order to facilitate this 
selection, the GEANT4 developers recommend to start cither from available examples 
of previously developed applications, or from the so-called predefined physics lists. 

Here we briefly describe the choice of models and main parameters used in our 
calculations. In MCHIT the energy loss of primary and secondary charged particles 
due to electromagnetic processes is described via a set of models called 'standard 
electromagnetic physics'. It accounts for energy loss and straggling due to interaction 
with atomic electrons as well as multiple Coulomb scattering on atomic nuclei. 

At each simulation step, the ionisation energy loss of a charged particle is 
calculated according to the Bethe-Bloch formula. The mean excitation potential of 
water molecules was set to 77 eV, i.e. to the value which better describes the set 
of available data on depth-dose distributions, both for proton and carbon-ion beams, 
in the range of therapeutic energies of 80 - 330 MeV per nucleon. For PMMA, the 
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mean excitation potential was set to 68.5 eV, a default value used in GEANT4 for 
this material. 

The bremsstrahlung processes for electrons and positrons were activated in the 
simulations along with the annihilation process for positrons. For photons, Compton 
scattering, the conversion into an electron-positron pair and the photoelectric effect 
were included. 

Two kinds of hadronic interactions are considered in the MCHIT model: (a) 
elastic scattering of hadrons on target protons and nuclei, which dominate at low 
projectile energies, and (b) inelastic nuclear reactions induced by fast hadrons and 
nuclei. The model of nucleon-nucleon elastic scattering is based on a parameterisation 
of experimental data in the energy range of 10-1200 MeV. At higher energies the 
hadron-nucleus elastic scattering is simulated within the Glauber model (GEANT4- 
Documents 2006). 

Overall probability of hadronic interactions for nucleons and nuclei propagating 
in the media is determined by the total inelastic cross section for nucleon-nucleus 
and nucleus- nucleus collisions. Parameterisations by Wellisch and Axen (1996) that 
best fit experimental data were used to describe the total reaction cross sections in 
nucleon-nucleus collisions. Systematics by Shen et al (1989) was used for the total 
nucleus-nucleus cross sections. 

For inelastic interactions of hadrons two groups of models are available in the 
GEANT4 toolkit (Agostinelli et al 2003): (a) the data-driven models, which are based 
on the parameterisations of measured cross sections for specific reaction channels, and 
(b) the theory-driven models, which are based on various theoretical approaches and 
implemented as Monte Carlo event generators. 

In the MCHIT model the inelastic interaction of low-energy (below 20 MeV) 
nucleons, including radiative neutron capture, were simulated by means of data 
driven models. Above 20 MeV the exciton-based precompound model was invoked 
(Agostinelli et al 2003). 

For hadrons and nuclei with energies above 80A MeV, we have employed the 
binary cascade model (Folger et al 2004). In this case exited nuclear remnants are 
created after the cascade stage of interaction. Therefore, appropriate models for the 
de-excitation process should be included in the simulation. The Weisskopf-Ewing 
model (Weisskopf and Ewing 1940) was used for the description of evaporation of 
nucleons from nuclei at excitation energies below 3 MeV per nucleon. The Statistical 
Multifragmcntation Model (SMM) by Bondorf et al (1995) was used at excitation 
energies above 3 MeV per nucleon to describe multi-fragment break-up of highly- 
excited residual nuclei. The SMM includes as its part the Fermi break-up model 
describing an explosive disintegration of highly-excited light nuclei. 

Various unstable nuclei produced in nuclear reactions were followed until their 
decay, in particular, the /3 + -decay leading to the emission of a positron. Within 
the GEANT4 toolkit such processes are simulated by using the tables of radioactive 
isotopes with the corresponding decay channels and their probabilities. 

4. Production of positron-emitting nuclei by protons and carbon ions 

In the case of protons propagating in PMMA the positron-emitting nuclei 11 C, 10 C 
and 15 are created mostly via (p,n) or (p,2n) reactions on 12 C and 16 O nuclei. In this 
case all the positron-emitting nuclei are fragments of target nuclei which are initially 
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at rest. In contrast, carbon beams produce n C and 10 C nuclei mostly via projectile 
fragmentation, while 15 nuclei are produced from the target 16 nuclei. 

The total /3 + -activity yields were measured for n C, 10 C and 15 by Parodi (2004) 
for protons and carbon ions stopped in PMMA phantoms. Since the half-lives of these 
nuclides are essentially different and well-known from the literature, one can build 
a mathematical model describing their decay and the total rate of /3 + -activity as a 
function of time. The measured time-dependence of the total /3 + -activity was fitted 
by Parodi (2004) with the yields of 11 C, 10 C and 15 considered as free parameters. 
The experimental yields of these nuclides obtained in such a way are listed in tables ^ 
and [21 along with their uncertainties. 

Calculational results of the MCHIT model for n C, 10 C and 15 are listed in 
tables H an d El for comparison. According to the MCHIT model, the yields of other 
/3 + -emitting nuclei with TWa > 1 min can be neglected. For example, the yields of 
13 N, 14 0, 17 F and 18 F nuclei together account for less than 5% of the total yield of 
/3 + -emitters produced by ~ 300A MeV 12 C beam in PMMA. The calculations were 
performed for the PMMA phantoms of the same size (300 mm x 90 mm x 90 mm) as 
used by Parodi (2004) in measurements. The beam profile in the transverse directions 
was assumed to be a Gaussian with the FWHM of 10 mm, while the beam energy 
spread was taken with the FWHM of 0.2%. 

The results of the FLUKA and POSGEN codes were reported by Parodi (2004) 
for several beam energies. Some of them are listed in tables ^ an( l HI respectively. 



Table 1. Calculated yields of positron-emitting nuclei (per beam particle, in %%) 
produced by 110, 140 and 175 MeV protons in PMMA phantom. Experimental 
data and FLUKA results by Parodi (2004) are shown for comparison. 





110 MeV 




140 MeV 




175 MeV 




MCHIT 


Experiment 


MCHIT 


Experiment 


FLUKA 


MCHIT 


Experiment 


n c 


1.83 


2.2 ±0.3 


2.64 


3.4 ±0.4 


2.67 


3.71 


4.7 ±0.7 


10 c 


0.11 


0.09 ±0.03 


0.20 


0.15 ± 0.03 


0.10 


0.31 


0.17 ±0.06 


16 


0.80 


0.80 ±0.15 


1.10 


1.23 ± 0.18 


1.23 


1.54 


1.6 ±0.3 



For protons stopped in the PMMA phantom, the yields of 15 O nuclei predicted by 
the MCHIT model show very good agreement both with the experiment and with the 
FLUKA results. Since n C and 15 are the most abundant and long-lived nuclides, an 
accurate description of their yields is crucial for the PET interpretation. As seen from 
table m the agreement between MCHIT results and the experiment is within ~ 20% 
accuracy. As compared with n C and 15 yields, the yield of 10 C is much lower for 
proton-induced reactions. The MCHIT results for 10 C are in good agreement with 
the experiment for 110 MeV protons, while for 140 and 175 MeV protons the theory 
overestimates 10 C yields as compared with the experimental results. 

As follows from table the predictions of the POSGEN model for n C and 15 
nuclides, reported by Parodi (2004) for 270. 55A MeV 12 C ions, are less accurate as 
compared with the MCHIT model in describing the experimental data at 259. 5 A 
MeV. Apparently, more extended benchmarking of the models against each other and 
experimental data is needed to make a firm conclusion about their performance. 

The yields of 10 C nuclei are overestimated by the MCHIT model also for the 
12 C beam, in particular, for the most energetic beam of 343. 46A MeV. The POSGEN 
predictions for the beam energy 270. 55A MeV are closer to the experimental results 
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Table 2. Calculated yields of positron-emitting nuclei (per beam particle, in %%) 
produced by 212.12A, 259. 5A and 343.46A MeV 12 C ions in PMMA phantom. 
Experimental data at the same beam energies and POSGEN results (calculations 
for 270.55A MeV 12 C) by Parodi (2004) are shown for comparison. 





212.12A MeV 




259.5A MeV 




343.46A MeV 




MCHIT 


Experiment 


MCHIT 


Experiment 


POSGEN 


MCHIT 


Experiment 


n c 


11.9 


10.5 ± 1.3 


16.83 


14.7 ± 1.6 


26.6 


25.25 


19.9 ±2.4 


io c 


1.97 


0.8 ±0.3 


2.79 


1.2 ±0.3 


1.96 


4.27 


1.5 ±0.3 


15 


2.38 


2.1 ±0.3 


3.69 


3.1 ±0.4 


10.0 


6.09 


5.0 ±0.4 



obtained at 259. 5A MeV. Nevertheless, we conclude that the MCHIT model is able 
to describe reasonably well the total yields of the most abundant positron-emitting 
nuclei, n C and 15 0, produced by proton and 12 C beams in PMMA. 

5. Spatial distribution of positron-emitting nuclei 

As shown above, the MCHIT model is quite successful in describing the absolute 
yields of /3 + -emitting nuclei. One can now study the spatial distributions of these 
nuclei. Such distributions were calculated for 110 MeV protons and 212. 12A MeV 
12 C ions together with the corresponding depth-dose distributions for these beams. 
The results are presented in figure ^ At these energies, protons and carbon ions have 
similar ranges in PMMA phantoms. In both cases the Bragg peaks are located at 
the depth of ~ 80 mm. As demonstrated in the previous section, the contribution 
from n C and 15 nuclei dominate for both the proton and carbon beams. However, 
the shapes and the absolute values of corresponding /3 + -activity distributions are very 
different. Protons produce 11 C, 10 C and 15 O fragments from the target carbon and 
oxygen nuclei. The cross sections for the 12 C(p,pn) 11 C and 16 0(p,pn) 15 reactions 
increase with decreasing proton energy up to 40-50 MeV and then rapidly fall down. 
For proton energies below 20 MeV there is no neutron emission at all. As shown in 
figure^ n C and 15 dominate, and their yields are slowly increasing as the projectile 
protons slow down. The production of 10 C is significantly suppressed because removal 
of two neutrons without destroying the remnant is quite improbable. As soon as the 
proton energy becomes below the neutron emission threshold, the production of /3 + - 
activity is ceased, so that no /3 + -emitters are created by protons within ~ 1 cm before 
the Bragg peak. 

In contrast to the proton irradiation, the maximum of /3 + -activity produced by 
carbon ions is located very close to the Bragg peak and thus clearly marks its position. 
The projectile fragments of interest, 11 C and 10 C, are created in peripheral nucleus- 
nucleus collisions, and have velocities which are very close to the projectile velocity. 
At same initial velocity, the ranges of energetic ions in matter R are proportional to 
A/Z 2 , where A and Z are the ion mass and charge. Therefore, if 11 C and 10 C are 
created at zero depth, their ranges are shorter as compared to the range of projectile 
ions: R( n C) ~ 11/12 x i?( 12 C) andii( 10 C) - 10/12 x R( 12 C). However, such 11 C and 
10 C nuclei can be produced at any depth within the range of primary nuclei, excluding 
only the last few millimetres of the beam range, where 12 C are not energetic enough 
for fragmentation. As discussed in details, in particular, by Fiedler et al (2006), the 
spread in production points leads to the spread of stopping points of fragments, as for 
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Figure 1. Calculated depth-distributions of deposited energy (dash-dotted 
curves) and /3+-activity (histograms) for (a) 110 MeV protons and (b) 212. 12A 
MeV 12 C nuclei in PMMA phantom. The distributions of n C, 10 C and ls O 
nuclei are shown by the long-dashed, dashed and dotted histograms, respectively, 
and their sum is shown by the solid-line histogram. The distribution of actual 
e~*~ -annihilation points is shown by the short-dashed histogram. 

example, for n C fragments this spread can be estimated as Ax - 1/12 x i?( 12 C). 

Despite this spread, all the projectile fragments are stopped before the 12 C Bragg 
peak. As shown in figure ^ the sharp fall-off in the /3 + -activity distribution clearly 
indicates the position of the Bragg peak. This is a very attractive feature of carbon-ion 
therapy which makes possible on-line monitoring of the ions stopping points. 

The distribution of 15 nuclei in 12 C induced reactions looks similar to that for 
the proton irradiations, since in both cases 15 O nuclei arc produced from the target 
16 nuclei. This is why the distribution is quite flat and falls down before the Bragg 
peak. Strictly speaking, some of 11 C, 10 C and 15 nuclei are produced off the target 
carbon and oxygen nuclei by secondary particles like p,n and a-particles created in 
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projectile fragmentation. The contribution of such processes is small, but still visible 
in figure Q] beyond the Bragg peak, as light secondary particles propagate further. 

From the above analysis we see that 11 C fragments are the most suitable nuclei 
for monitoring the energy deposition in carbon-ion therapy. They have the largest 
yield and longest life time (Ti/ 2 — 20.39 min) as compared with other /3 + -emitters. 
Thus, mostly 11 C nuclei will survive for 10-20 min after stopping, and then emit the 
positrons. 

6. Distribution of positron annihilation points as measured by a PET 
scanner 

As explained in section [21 one can measure the depth distribution of /3 + -activity via 
Positron Emission Tomography (PET). Our present study is focused on the ability of 
the GEANT4 toolkit to describe the production of positron-emitting nuclei rather 
then on modelling various aspects of PET measurements. As has been already 
demonstrated by Jan et al (2004), Schmidtlein et al (2006) and Lamare et al (2006), 
the GEANT4 toolkit is very successful in Monte Carlo simulations of commercial PET 
scanners. 

Here we consider only one aspect of the PET monitoring method, i.e. how model 
predictions are affected by a finite resolution of PET scanners. In our model the PET 
signal is generated in the following simplified way. First, the decays of all unstable 
nuclei are simulated at their stopping points. Second, all the emitted positrons are 
traced up to their annihilation, and, finally, the distribution of e + e~-annihilation 
points is obtained. It is instructive to note that positrons with energy of a few 100 
keV have quite a long mean free path with respect to annihilation in materials like 
water, namely \ ann ~ 10 cm (a ann ~ 0.5 b) at normal water density. Therefore, the 
annihilation happens only when the positrons have already slowed down. 

Third, these distributions are convoluted with Gaussian spreading functions in 
order to mimic a realistic response of a PET imaging system. We choose the FWHM 
in the range of 8-^10 mm. Such spatial resolution was reported by Ponisch et al (2003) 
for their 3D PET reconstruction algorithm. Results of our simulations for proton and 
carbon- ion beams are shown separately in figures [21 and [3] 

Since the measured /3 + -activity distribution depends on the spatial resolution of 
the PET scanner, this simple procedure makes possible to compare the calculated 
/3 + -activity profiles with those obtained from experiment, as shown in figures [21 and [21 
The distributions of e + annihilation points generated with the MCHIT model are 
shown in figure [21 for 62 and 110 MeV protons in PMMA phantoms. In both cases 
the distributions are very wide. When the finite PET resolution is taken into account 
via the convolution of these distributions with a Gaussian of FWHM = 8 mm, only 
small changes in the shape of the distributions are obtained. Nevertheless, this leads 
to better description of experimental data by Oelfke et al (1996). The agreement 
with the experiment is very good for 62 MeV protons. On the other hand, some 
disagreement with the experimental data is found for 110 MeV protons, i.e. the 
measured distribution is more flat as compared to the calculated one. In figure the 
MCHIT results are compared with the PET data for carbon ions in PMMA phantoms 
obtained by Ponisch et al (2004) for 212. 12A MeV, by Parodi (2004) for 259. 5A MeV 
and by Pawelke et al (1997) for 292. 5A MeV. Some monitoring devices were placed 
in front of the PMMA phantoms during measurements while the MCHIT calculations 
were performed for pure PMMA phantoms without any additional elements. To 
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Figure 2. Calculated depth-distributions of deposited energy (dash-dotted 
curve) and positron annihilation points (histograms) for (a) 62 MeV and (b) 
100 MeV protons in PMMA phantom. The distribution of actual e+ annihilation 
points is shown by the dashed histogram, while the distribution which accounts 
for a finite spatial PET resolution of FWHM = 8 mm is shown by the solid 
histogram. Points show experimental data (Oelfke et al 1996). 



correct for this difference, the experimental data in figure |3 were shifted by 3 - 6 
mm to ensure the same position of the Bragg peak in PMMA as reported by Parodi 
(2004). 

For carbon beams, the consideration of the actual PET resolution is crucial for 
proper description of the experimental data. In this case the actual distribution of e + 
annihilation points consists of a sharp peak on a flat plateau. In order to investigate 
the effect of the PET scanner resolution, the calculated distributions were convoluted 
with the Gaussians of FWHM = 8 and 10 mm, resulting in the distributions shown 
in figure [3J The calculation with FWHM = 8 mm shows a good agreement with 
the experimental data at 212. 12A MeV, while for 259. 5A MeV a better agreement is 
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Figure 3. Calculated depth-distributions of deposited energy (dash-dotted 
curve) and positron annihilation points (histograms) for 12 C beam with energy 
212A MeV (a), 259.5A MeV (b) and 292.5A MeV (c) in PMMA phantom. 
The distribution of actual e + annihilation points is shown by the long-dashed 
histogram, while the distributions which account for a finite spatial PET 
resolution of FWHM = 8 and 10 mm are shown by the solid and short-dashed 
histograms, respectively. The experimental data are shown by light solid lines: 
(a) ( Ponisch et al 2004) (b) (Parodi 2004) and (c) (Pawelke et al 1997). 

achieved for FWHM = 10 mm. The MCHIT calculations at 292. 5A MeV produce a 
lower plateau as compared with experimental data by Pawelke et al (1997). Since the 
measurements in this case were made promptly after irradiations, the contribution 
of short-lived isotopes populating the plateau should be bigger compared to the 
measurements at 212. 12A and 259. 5A MeV. One can note that the calculations by 
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Pawelke et al (1997) also show the peak-to-plateau ratio, which is higher than found 
in their measurements. 

7. Comparison of dose and [3 + activity distributions in different materials 

In order to use the PET method for dose monitoring, the employed calculational 
tool must well describe both the dose and /? + -activity distributions. As the In- 
activity distributions have been already considered above, we shell study now the 
corresponding dose distributions calculated with MCHIT, which one can then verify 
with corresponding experimental data. 

The data on dose distributions in PMMA phantoms are rather scarce. Matsufuji 
et al (2003) have performed such measurements for beams of 4 He, 12 C, 20 Ne, 28 Si and 
40 Ar nuclei. In figure 0] the calculated dose distribution for 279.2 A MeV 12 C nuclei 
in PMMA is compared with the experimental data of Matsufuji et al (2003). 

As known, the charged particle range in material of a given chemical composition 
is inversely proportional to the density of electrons, and hence, to the matter density p. 
In particular, it is quite common to present the depth-dose distributions as functions 
of the areal density x ■ p. Then, if dE /Ax is known at one density, it can be easily 
calculated for other density by rescaling the depth x. Strictly speaking, this scaling 
procedure is fully justified only for a given material at various densities. However, one 
can try to use it also for different tissue-like materials. 

In figure 0] the scaling property was used to compare the calculated dose and 
/3 + -activity distributions in water with those in PMMA and in dense bone tissue. The 
elemental composition of the bone tissue was taken in the following mass fractions: H 
- 6.4%, C - 27.8%, N - 2.7%, O - 41%, Mg - 0.2%, P - 7%, S - 0.2%, Ca - 14.7%. 

In order to get the same position of the Bragg peak, the depth in PMMA and 
bone were stretched by factors of 1.16 and 1.74, respectively, which are close to the 
density ratios, ppMMA/pwatcr = 1.18 and pbone/pwatcr = 1-85. In order to account 
for the changed bin size, the histograms were rescaled by the inverse factors, 1/1.16 
and 1/1.74, respectively. Plotted as functions of water- equivalent depth, the dose 
and /3 + -activity profiles in water and PMMA are in very good agreement with each 
other despite of the difference in the chemical compositions of PMMA (C5H8O2) and 
water (H 2 0). However, the calculated position of the Bragg peak is ~ 4 mm deeper 
in water compared to the depth-dose profile measured by Matsufuji et al (2003). As 
discussed by Gudowska et al (2004), the position of the Bragg peak in experiment 
was determined with an accuracy of ±3 mm. This discrepancy can be attributed to 
the uncertainty in the initial energy of the carbon beam. In order to reach agreement 
with the experimental distributions the data by Matsufuji et al (2003) were shifted by 
4.2 mm towards the larger depth, see figure This is exactly the same shift, which 
was used by Gudowska et al (2004) to obtain the agreement between the data and the 
results of their SHIELD-HIT code. 

As shown in figure 01 the /3 + -activity distributions in water and PMMA can be 
successfully transformed from one material to another by rescaling. However, this 
procedure fails for another pair of materials, water and bone tissue. As one can see 
in figure the scaling transformation does not correctly reproduce the Monte Carlo 
simulation results for the energy deposition in bone tissue. The-peak-to-plateau ratio 
for the /3 + -activity distribution in bone tissue is lower compared to water due to the 
difference in their elemental compositions. As explained above, this peak is mainly 
generated by the 11 C nuclei produced from the projectile 12 C nuclei. These nuclei 



Positron emitting nuclei in proton and heavy-ion therapy 



13 




20 40 60 80 100 120 140 160 180 200 

water-equivalent depth (mm) 



Figure 4. (a) Dose distributions as functions of water-equivalent depth 
calculated for 279. 2A MeV 12 C nuclei in water (solid-line histogram), PMMA 
(dashed histogram) and bone (dotted histogram). Experimental data for depth- 
dose distributions in PMMA (Matsufuji et al 2003) are shown by points, (b) 
Distributions of positron-emitting nuclei produced in these materials. 

are produced less frequently in collisions with target nuclei of bone tissue, which are, 
on average, more heavy compared with water. In addition, a larger part of the total 
/3+-activity (~ 10%) in bone tissue is due to 13 N, 14 0, 17 F, 18 F and 30 P nuclei, which 
are produced in target fragmentation reactions. This contribution, which also includes 
the reactions induced by secondary nucleons beyond the Bragg peak, is characterised 
by a flat depth-distribution. 

We believe that after appropriate validation, our Monte Carlo approach can 
be used for calculating /3 + -activity profiles in various materials, as an alternative 
method to the less accurate scaling transformation. This approach should be especially 
powerful in the case of strongly inhomogeneous medium as human body where soft 
tissues are intermittent with voids and bones. 
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8. Conclusion 

In the present work the validation of the GEANT4 toolkit for calculations in particle 
therapy is extended to the /3 + -activity distributions induced by proton and 12 C 
beams. Despite the fact that these /3 + -activity distributions differ in shape from 
the corresponding dose distributions, there exists a strong correlation between them, 
which can be used for the dose and range monitoring by means of the PET scanners. 
Since human tissues and organs have various nuclear compositions, only detailed 
Monte Carlo simulations of proton and carbon-ion therapy are able to account for 
the differences in /3 + -activity distributions for different kinds of tissues, e.g. soft 
tissues and bones 

Our MCHIT model based on the GEANT4 toolkit is able to predict with 
reasonable accuracy (1) the total yields of n C and 15 nuclei, which are the most 
abundant /3 + -emitters produced by proton and 12 C beams, (2) distributions of e + 
annihilation points as measured by PET scanners with a realistic spatial resolution, 
(3) the deviations from simple density scaling in the dose and /3 + -activity distributions. 
By this study, the ability of the GEANT4 toolkit to correctly simulate the dose 
and /3 + -activity distributions in the proton and carbon ion therapy has been clearly 
demonstrated. Therefore, we suggest to use the MCHIT model for realistic calculations 
of /3 + -activity profiles and dose distributions in proton and carbon-ion therapy. 
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